def get_mesh_info(outputdir="mesh_output"):

    import string
    import numpy as np

    # the three files we need
    parameter_file_name = outputdir + "/mesh_params.dat";
    parameter_file      = open(parameter_file_name,'r');
    node_xyz_file_name  = outputdir + "/mesh_xyz_node.dat";
    node_xyz_file       = open(node_xyz_file_name,'r');
    tnode_file_name     = outputdir + "/mesh_tnode.dat";
    tnode_file          = open(tnode_file_name,'r');

    # grab number of elements and number of nodes
    linestring = parameter_file.readline(); linelist = linestring.split();
    NumElems = int(linelist[0]);
    linestring = parameter_file.readline(); linelist = linestring.split();
    linestring = parameter_file.readline(); linelist = linestring.split();
    linestring = parameter_file.readline(); linelist = linestring.split();
    NumNodes = int(linelist[0]);
    parameter_file.close();

    # read-in all nodes in xyz format
    node = np.zeros([NumNodes,3]);
    for i in range(0,NumNodes):
        linestring = node_xyz_file.readline();
        linelist = linestring.split();
        node[i,0] = float(linelist[0]);
        node[i,1] = float(linelist[1]);
        node[i,2] = float(linelist[2]);

    # read-in all triangles
    tnode = np.zeros([NumElems,3],dtype=int);
    for i in range(0,NumElems):
        linestring = tnode_file.readline();
        linelist = linestring.split();
        tnode[i,0] = int(linelist[0])-1;
        tnode[i,1] = int(linelist[1])-1;
        tnode[i,2] = int(linelist[2])-1;

    return NumNodes,NumElems,node,tnode

#--------------------------------------------------------------------------
if __name__ == "__main__":

    # run mesh generator
    import os;
    os.system('mesh.exe &> /dev/null');

    # grab mesh information
    NumNodes,NumElems,node,tnode = get_mesh_info();

    # grab quadrature points
    import gauss_quad_sphere as gqs;
    quad_pts,quad_wgts,NumQuadPts = gqs.get_gauss_quad_sphere(6);

    # compute surface area of each element
    import numpy as np;
    onethird = 1.0/3.0;
    surf_area = np.zeros(NumElems);
    total_area = 0.0;
    for i in range(0,NumElems):
        t = tnode[i,0:];
        x1 = node[t[0],0]; y1 = node[t[0],1]; z1 = node[t[0],2];
        x2 = node[t[1],0]; y2 = node[t[1],1]; z2 = node[t[1],2];
        x3 = node[t[2],0]; y3 = node[t[2],1]; z3 = node[t[2],2];

        a = -y2*z1 + y3*z1 + y1*z2 - y3*z2 - y1*z3 + y2*z3;
        b =  x2*z1 - x3*z1 - x1*z2 + x3*z2 + x1*z3 - x2*z3;
        c = -x2*y1 + x3*y1 + x1*y2 - x3*y2 - x1*y3 + x2*y3;

        xcenter = onethird*(x1+x2+x3);
        ycenter = onethird*(y1+y2+y3);
        zcenter = onethird*(z1+z2+z3);

        tmp = 0.0;
        for k in range(0,NumQuadPts):
            xi  = quad_pts[k,0];
            eta = quad_pts[k,1];
            omega = quad_wgts[k];
            xx  = xcenter + xi*(x2-x1) + eta*(x3-x1);
            yy  = ycenter + xi*(y2-y1) + eta*(y3-y1);
            zz  = zcenter + xi*(z2-z1) + eta*(z3-z1);
            r3  = (np.sqrt(xx**2+yy**2+zz**2))**3;
            tmp = tmp + omega*(a*xx + b*yy + c*zz)/r3;

        surf_area[i] = tmp;
        total_area = total_area + tmp;

    print(" ")
    for k in range(0,NumElems):
        print(surf_area[k])
    print(" ")
    print(np.abs(4.0*np.pi-total_area)/(4.0*np.pi))
    print(np.min(surf_area))
    print(np.max(surf_area))

#--------------------------------------------------------------------------
